%% SOLVES HJB IN LIMITED-COMMITMENT CONTRACT
function V = limited_hjb_solver(V00, Beta, f, mu_v, vol_phi, params, grid_params, num_params)
%% PARAMTERS
temp = num2cell(params);
[rho, r, ra, alpha, sigma, lambda, M, eta_min, eta_max] = temp{:};

[n_grid, ~, PHI_mat, ~, dphi] = grid_params{:};

[max_iter, D, thr] = num_params{:};


%% v-term
mat_v = mu_v;

AA_v = spdiags(mat_v, 0, n_grid,n_grid);

%% AA matrix for v-term and first derivatives
AA = AA_v;

%% Second derivatives
vol2_phi = vol_phi.^2/(2*dphi^2);

up_mat_vol_phi = [0;vol2_phi];

cent_mat_vol_phi = -vol2_phi*2;

low_mat_vol_phi = [vol2_phi(2:end);0];

BB_phiphi =spdiags(cent_mat_vol_phi,0,n_grid,n_grid)+...
    spdiags(low_mat_vol_phi,-1,n_grid,n_grid)+...
    spdiags(up_mat_vol_phi,1,n_grid,n_grid);

BB = BB_phiphi;

%% "Cross derivatives" with Scaling Variable
cross_phi = vol_phi.*Beta/(dphi);
up_cross_phi = max(0,cross_phi);
cent_cross_phi = - max(0,cross_phi) - max(0,-cross_phi);
low_cross_phi = max(0,-cross_phi);

up_mat_cross_phi = [0;up_cross_phi];
cent_mat_cross_phi = cent_cross_phi;
low_mat_cross_phi = [low_cross_phi(2:end);0];

CC_phi = spdiags(up_mat_cross_phi,1,n_grid,n_grid)+...
         spdiags(cent_mat_cross_phi, 0, n_grid, n_grid)+...
         spdiags(low_mat_cross_phi,-1,n_grid,n_grid);    

CC = CC_phi;

%% New value function
P = AA + BB + CC;

Q = ((r + 1/D)*speye(n_grid) - P);

V = Q\(f + V00/D);
